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Abstract 

We consider a process in which there are two types of particles, A and on an 
infinite one-dimensional lattice. The particles hop to their adjacent sites, like the 
totally asymmetric exclusion process (ASEP), and have also the following interactions: 
A + B B + B and B + A B + B, all occur with equal rate. We study this 
process by imposing four boundary conditions on ASEP master equation. It is shown 
that this model is integrable, in the sense that its A^-particle S-matrix is factorized 
into a product of two-particle S-matrices and, more importantly, the two-particle S- 
matrix satisfy quantum Yang-Baxter equation. Using coordinate Bethe-ansatz, the A^- 
particle wavefunctions and the two-particle conditional probabilities are found exactly. 

Further, by imposing four reasonable physical conditions on two-species diffusion- 
reaction processes (where the most important ones are the equality of the reaction rates 
and the conservation of the number of particles in each reaction), we show that among 
the 4096 types of the interactions which have these properties and can be modeled by 
a master equation and an appropriate set of boundary conditions, there are only 28 
independent interactions which are integrable. We find all these interactions and also 
their corresponding wave functions. Some of these may be new solutions of quantum 
Yang-Baxter equation. 
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1 Introduction 



For non-equilibrium systems in low dimensions, an understanding can often be gained by 
studying rather simple models []T1-0|. One of the important examples of these systems are 
reaction-diffusion processes on a one-dimensional lattice, which their dynamics are fully 
specified by their master equation 0-0. In some cases, it is possible to solve the master 
equation exactly. In recent years, there has been enormous progress in the field of exactly 
solvable non-equilibrium processes. These developments were mainly triggered by the ob- 
servation that the Liouville operator of certain (l-l-l)-dimensional reaction-diffusion models 
may be related to Hamiltonians of previously known quantum spin systems 0, |0. 

One of the simplest examples of reaction-diffusion processes are Asymmetric Simple 
Exclusion Processes (ASEP) 0, 0, which has been used to describe various problems 



in different fields of interest, such as the kinetics of bipolymerization[|r^, dynamical models 
of interface growth [|1^, and traffic models[|l^]. The totally ASEP model has been solved 
exactly by imposing the appropriate boundary condition on the probabilities appear in the 
master equation ||14||. The totally ASEP model describes a process in which each lattice site 
can be occupied by at most one particle and the particle hops with rate one to its right 
neighboring site if it is not already occupied, otherwise the attempted move is rejected. 

There are some other interesting and more complicated processes which can be solved 
by the method developed in [|14|, namely by choosing a suitable boundary condition for the 
master equation. For example, it has been shown that the so called "generalized totally 
ASEP model" can be solved exactly by this method |]T5[. In this model, even if the right 
neighboring site of a particle is occupied, the particle hops to the next right site by pushing 
all the neighboring particles to their next right sites, with a rate which depends on the 
number of right neighboring particles. This model has been further generalized in [|T^ by 
considering both the right and left hopping of the particles. 

In this paper we are going to consider a class of integrable models in which there are two 
species of particles which can hop to their right neighboring sites if those are not occupied, 
and also the particles interact with each other if they are in adjacent sites. The details of 
this nearest-neighboring interaction depends on the specific considered model (see ||T7[]-[p!S 



for some recent works in two- and three-species reaction-diffusion processes). The impor- 
tant point in integrable reaction-diffusion processes with more than one type of particle is 
that, as we will show, the two-particle S-matrix of the reaction, which specify the A^-point 
functions, must satisfy the Quantum Yang-Baxter Equation (QYBE). Therefore, as we ex- 
pect, the number of integrable models, in the sense that its A^-particle S-matrix can be 
factorized into a product of two-particle S-matrices, is very few. In this paper we will find 
all two-species integrable reaction-diffusion processes which have the following properties: 
1. the particles hop to their right neighboring sites if these sites are not occupied. 
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2. the interaction occurs only between nearest-neighbor particles, 

3. the particles can be annihilated or created, with the only restriction that the total number 
of particles is fixed, 

4. all the interactions, including diffusions, occur with the same rate. 

We show that among the 2^^ = 4096 types of the interactions which have the above- 
mentioned properties and can be modeled by a master equation and an appropriate boundary 
condition, there are only 42 interactions which are integrable (their two-particle S-matrices 
satisfy the QYBE), and from these 42 interactions, only 28 of them are independent. Some 
of these may be new solutions of QYBE. 

The plan of the paper is as following. In section 2, we introduce the first kind of this 
interactions, which was our initial motivation in this work, in which besides the usual hop- 
ping, the two types of particles interact as: A + B ^ B + B and B + A B + B. We 
show that this interaction can be modeled by the usual master equation of ASEP and four 
boundary conditions. We also show that the model is integrable. Note that one can look 
at this model (see eq.(P) as a simple one-dimensional model of spread of disease. If we 
consider A particles to be the healthy individuals and the B particles the diseased ones, 
then we expect that when A and B particle are near to each other, healthy one may become 
diseased (in other words B transmits disease to A). In section 3, we compute the exact two- 
particle conditional probabilities of this interaction and study the long-time behavior of this 
probabilities. And finally in section 4, we investigate the class of integrable models which 
have the four above-mentioned properties and deduce that there are 28 different models, 
which the totally ASEP model and our first model introduced in section 2, are two of them. 

2 AB BB and BA BB reaction diffusion process 
2.1 The master equation 

Suppose there are N particles of two types A and B on an infinite one dimensional lattice, 
with interactions 

AB BB, ^ ' 

BA BB, 

all occur with equal rate, which can be scaled to one. In eq.([l|), we denote the vacancy by 
notation 0. The basic quantities we are interested in are the probabilities Poj^aj^.-.^aj^ (xi, X2, ■ ■ 
■,X]\f;t) for finding at time t a particle of type ai at site Xi, a particle of type 02 at site X2, 
etc. Each a, can be A or B. Following |]T^, we take these functions to define probabilities 
only in the physical region xi < X2 < ■ ■ ■ < xn- The surfaces where any of the two adjacent 
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coordinates are equal, are the boundaries of the physical region. In the subset of the physical 
region where Xj+i — Xj > 1, Vi, we have only hopping of the particles and therefore the master 
equation can be written as 

+-Pai,-,aiv(a;i,a;2, ■ ■ -^xn - (2) 

As is clear from eq.(0), when Xj+i = Xj + 1 for some i's, the one or more of the probability 
functions go out from the physical region and therefore the eq.(0) has to be supplemented 
by some boundary conditions. The particular choice of the boundary condition depends on 
the details of the interactions of particles. For example, it can be shown that in the totally 
ASEP model, the suitable boundary condition is |]14| , 



P{x,x) = P{x,x + 1), Wx, (3) 

in which the time variable and also all the other coordinates have been suppressed for 
simplicity. The master equation and the boundary condition (|^) replace the very large 
number of equations which one should write by considering the multitude of cases which 
arises in different possible configurations. 

To model the interaction (|l]), we claim that the suitable boundary conditions are : 

PAAix, x) = Paa{x, x + l), (4) 

Pbb{x, x) = Pbb{x, x + 1) + Pab{x, x + 1) + Pba{x, x + 1), (5) 

Pab{x, x) = Pba{x, x) = 0, (6) 

where we have again suppressed the positions of all the other particles. By looking at (|l]), 
it is obvious that if we have only A particles , the process is exactly the same as totally 
ASEP. It is the reason of appearing eq.(|D which is the same as (|]). To justify the other 
three boundary conditions, we provide a few examples in the two- and three-particle sectors, 
instead of giving a general proof. 

First we consider the two-particle sector, for example Pba{x, x + 1). By master equation 
(H) we have 

d 

^Pba{x, x + l)= Pba{x - 1, x) + Pba{x, x) - 2Pba{x, x + 1). (7) 
ot 

Using d^), (0) reduces to 
d 

g^PBAix, x + l) = Pba{x - 1, x) - 2Pba{x, x + 1). (8) 
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This is exactly what we expect, as the source of configuration (05^40) is (-60^40) and its 
sinks are two configurations (050^4) and (05S0). The second example is Pbb{x,x + 1). 
Using again the master equation and the boundary condition we obtain 

d 

— Pbb(x, a; + 1) = Pbb{x - 1, x) + Pab{x, a; + 1) + Pba{x, a; + 1) - Pbb{x, x + 1). (9) 
at 

This equation also predicts the true source and sink terms, because {^BB^) has three sources 
{B(l)B(!i), (0A50), and (0-BA0) and only one sink (0505). As a three-particle sector example, 
let us consider the most nontrivial case Pbbb{x,x + l,x + 2). Using (|^) and (H), we find 

§iPbbb{x, X + 1,x + 2) = Pbbb{x - 1, a; + 1, a; + 2) + Pbbb{x, x,x + 2) + 
Pbbb{x,x + 1,x + 1) - 3Pbbb{x,x + l,a; + 2) 
= Pbbb{x -1,x+1,x + 2) + 2Pbab{x, x + 1,x + 2) + 
Pbba{x, a; + 1, a; + 2) + Pabb{x, a; + 1, a; + 2) - Pbbb{x, x + 1,x + 2). 

This is also the true equation, because the configuration {^BBB^) has five sources namely 
{BiDBB^),2{^BAB^), {%BBA%), and (0A550) and one sink {%BB%B). The reason of ap- 
pearing the factor 2 in (05AB0) is that {^BA ■ ■■) can go to (055 ■ ■■) and also (■ ■ -^50) 
can go to (■ ■ -550). 



2.2 The Bethe ansatz solution 

We now try to solve the master equation @) with boundary conditions ( ^-(^ by Bethe 
ansatz method. If we define ai,---,aNi.Xi, ■ ■ -yXj^f) through 

Pai, -,aN{xi, ■ ■ xn; t) = e"'^*^„i,...,„j^(xi, ■ ■ -, Xn), (11) 

and substitute (|lT]) in master equation (||) and boundary conditions (|D-(|D, we find 



and 



^ai,-,aiv(2^1 -l,X2,- ■ ■,Xn) -\ h ,...,ajv (^^l ' ^2, ■ ■ ■,Xn -I) 

= {N - eN)'^aj_,-,aN{Xi,X2,- ■ ■,Xn), (12) 

'^Aa{x,x) = "^AAix^X + 1), 

^BBix,x) = ^BBix,X+l) + '^ABix,X + l) + <ilBAix,X+l), (13) 
'^Ab{x,x) = ^BAix,x)=0. 



Following pO[, it becomes easier if we use a compact notation as follows: ^' is a tensor of 
rank N with components '^ai, - -,aNixi, ' ' ',Xn)- Then the boundary conditions (|l^) can be 
written as 

^(•••,e,e,---)=bfc,fc+i*(---,e,e+l,---), (14) 
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(15) 



where h^^k+i is the embedding of b (the matrix derived from ([T3|)) in the location k and 
k + 1: 

bfc.fe+i = 1 ® • ■ -(8) ^b^ ® ■ ■ • ®1, 
k,k + l 

with 

/ 1 \ 



Vo 1 1 i; 

To solve eq.(|12|), we write the coordinate Bethe ansatz for ^ in the form: 



(16) 



(17) 



where x and p denote the A^-tuples coordinates and momenta, respectively, the summation 
runs over all the elements of permutation group, and A^-'s (tensors of rank A^) are coefficients 
that must be determined from boundary condition (0). Inserting (|T^ into ([T^), yields: 

J2 A,e*"(P)-"(e-^"(P^) + ■ ■ ■ + e-'^^P''^ + - N) = 0, (18) 



from which one can find the eigenvalues eN as: 



AT 
k=l 

To find the coefficients A^., we insert the wavefunction ( |I7D in ([I^) , which yields 

We note that the first part of the above equation is symmetric with respect to interchange 
of pk and Pk+i, so if we symmetrize the bracket with respect to this interchange, we obtain 

(1 - e^'^(^''=+^)b,,fc+i)A^ + (1 - e^'^(P^)b,,fc+i)A^,, = 0, (21) 

where is an element of permutation group which only interchange and Pk+i, and aa^ 
stands for the product of two group elements a and cxfc. Thus we obtain: 

^aa^ = Sk,k+lior{pk), (T{pk+l))A^, (22) 

where 

Sfc,fc+l(2;i,22) = 1 ® ■ ■ S{Zi,Z2) ■ • ■ 01, 

^ ^ ' (23) 

k,k + l 



■'bfc ,fc + l)Acr 



0. 



(20) 
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and 

S{z,,Z2) = -{l-z^h)-\l-z,h). (24) 

In the above equations, we have denoted e'^'^'' by z^. In this way, we can calculate all the 
coefficients A^-'s in term of Ai by using eq.(^), and it seems that the problem is solved for 
arbitrary boundary condition b matrix. But this is not the case and we should moreover 
check the consistency of the solutions, which is highly nontrivial and depends on the details 
of the interaction, i.e. the elements of the b matrix . 



To see this, let us find two coefficients Ao-^ 0-20-1 and Ao-2o-io-2- Note that aia2cria.nd 02^102 
are equal as the elements of permutation group, that is both of them when act on (pi, P25 Ps, P45 
•) will result in 

(Pl,P2,P3,P4, •••) (P3,P2,Pl,P4, •••)• 

So we should have 



A 



(Tl(T20"l 



A 



(J2<Tl(T2 • 



(25) 
(26) 



Now using (0), we have: 
A 



^ia2<T3 = 5i2(aia2(pi), (Ti (12(^2)) A0102 = '5l2(p2,P3)A^^^2 

= 5'i2(p2,P3)5'23(cri(p2),0-l(P3))A^i = 5*12 (^2, P3) 5*23 (Pl , P3) )Ao 
= 5'i2(p2,P3)5'23(Pl,P3)5'l2(pi,P2)Ai, 



(27) 



and in the same way 



A020-10-2 — 

5'23 (Pl , P2 ) 5*12 (Pl , P3 ) 5'23 (P2 , P3 ) A 1 . 



Therefore, (^) yields: 

5'l2(P2,P3)5'23(Pl,P3)5'l2(Pl,P2) = 5'23 (pi , ^2) 5'l2 (Pl , P3) 5'23 (P2 , P3) , 



(28) 



(29) 



which is the familiar Quantum Yang-Baxter equation. Therefore , we must check if the 
S-matrices defined in (|^) and (p^ , with b from eq. ([T6D , satisfy (^) or not. 
Using (|T6]) , it can be shown that (0) is equal to 



S{z^ w) 



1 



1 — w 






1-z 
1-z 









(30) 



z — w z — w 1 — w 

and if we define z = e^^^ ,w = e*^^, and t = 6**^^, the eq. (^9]) can be written as 

{S{w, t) ® 1)(1 ® S{z, t)){S{z, w) ® 1) = (1 ® S{z, w)){S{z, t) S{w, t)). (31) 

Now it is not too difficult to show that the matrix ( pO|) really satisfy eq.(^ll), and therefore 
the solutions of interactions (|T]) are wavefunctions 
found by eq. (^2]) . 



with the coefficients that can be 
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3 The two— particle conditional probabilities 



In this section we want to calculate the two-particle conditional probabilities 
P{ai, a2, Xi, X2; t\i3i, P2, Vi, 112] 0), which is the probability of finding ai at site Xi and particle 
a2 at site X2 at time t, if at time t = we have the particle Pi at site yi and particle P2 at 
site y2- As has been discussed in and , we should take a linear combination of eigen- 
functions P{xi, X2) (from (|TT]) and ([T^)) with suitable coefficients, to find these two-particle 
conditional probabilities. Therefore, 



(x;t|/3,y;0) 



/ Paa \ 

Pab 
Pba 
V Pbb j 
I f{pi,P2)e~'^^'^ixi,X2)dpidp2 



(32) 



(2^) 



^ f e-^2tg-ip.y 





( ^ 




( ^ 


> 




< 


b 


e^iP^^^+p^^^) + Si2ipi,P2) 


b 




> dpidp2 




c 




c 








[d ) 




[d ) 







In the above expansion Pa^„2(x; t|/5,y; 0) is P(ai, 0^2, Xi,X2; t\Pi,P2, 2/1,1/2; 0), f{pi,P2) is the 
coefficient of expansion which in the second equality we choose it (2^^"*^'^' ^2 = 



-ipi 



(see (0)). ^(xi,X2) is the two-particle wave function where from ( p!7D and 



(|2^ ) we obtain 



(33) 



In the two-particle sector, Ai is a 4-column vector whose components must be determined 



by initial conditions, and 5*12(^1,^2) can be read from (|30D : 



'^12(^1,^2) 



/ s' \ 
0-100 
0-10 
\ s s s' ) 



(34) 



with 



oW2 



(35) 
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Let us first calculate eq.(0) irrespective of initial conditions , that is for arbitrary a, b, c, d 
Substituting (H) in (H), we find: 

/ a{F,{t) + F,{t)) \ 

biF.it) - F,it)) 
c(Fi(t)-F2(t)) 
V d{F^{t) + F,{t)) + {b + c)F^{t) J 



( Paa \ 
Pab 
Pba 
V Pbb ) 



(x;t|/?,y;0) 



(36) 



in which 



Fi{t) 
F2{t) 
Fsit) 
F,it) 



1 



(27r)2 
1 

1 
1 

(27r)2 



g-':2tg'«P-(x-y) 



dpidp2, 



e»pi - 1 
1 - e'P^ , 



.ei(p-^~p-y)dp^dp2 



(37) 
(38) 
(39) 
(40) 



e*pi - 1 

where in the above equations we have suppressed the x and p dependence of Fj's for simplicity 
and p = {p2,Pi)- To avoid the singularity in s and s', we set pi ^ pi + ie , and then by some 
simple contour integration we find 



i^2(0) = F3(0) = F4(0) = 0, 



(41) 



and at t 7^ 0, 



Fi{t) 
F2{t) 
Fsit) 

F,{t) 



-2t f^i-yi t="2-;/2 
i^i-yi)- (2^2-^2)! ' 

-2t P2-yi pi-a2 



(2:2-1/1)! (xj-y2)!' 



t^l-!/2 



(42) 



2t ( t^l-I/2 + l _ t^l-!/2 \ 

\{xi-y2+iy. (2^1-^2)!/ ^-^0 



P2-J/1 + ''' 

{x2-yi+ky. ■ 



Now we can study the different initial conditions. 



a) Case of /3i = /32 = A 



If at t = 0, both particles are A type, then our initial condition is 



/ Paa \ 
Pab 
Pba 
V Pbb J 



(x;0|AAy;0) 



/ ^xi,yi^X2,y2 ^ 







(43) 
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Using (|36|) and we find 

a=l,6 = c = rf = 0, (44) 

and therefore 

P^a(x; t\A, A, y; 0) = F,{t) + F^{t), (45) 
and all other P's are zero. Note that eq.(^5|) is exactly the same conditional probability that 



has been found in for simple ASEP model, as we expect, 
b) Case of l3i = A,l32 = B 

In this case the only non-zero element of conditional probability, at t = 0, is Pab = 

b = l,a = c = d = 0, (46) 



^xi,yiSx2,y2- Therefore one finds 



and therefore 

PAB{^;t\A,B,y;0) = Fi(t) - F2(t), 

PBB{^;t\A,B,y;0) = F^it), (47) 

and Paa = Pea = 0, which is consistent with our processes (|l]). It can be also checked that 
the conservation of probability holds, 

oo 2:2 — 1 

E iPAB + PBB){^;t\A,B,y;0) = l, (48) 

X2=y2 xi=yi 

for arbitrary yi,y2 and t. 

c) Case of Pi = 3,(52 = A 

In this case, the final result is : 

PBA{^;t\B, A, y;0) = F^{t) - F2{t), 

PBB{^;t\B,A,y;0) = F,{t) (49) 

and Paa = Pab = 0. 

d) Case of /3i = P2 = B 

In this case the final result is the same as the case 3.1, as we expect, 

Pbb(x; t\B, B, y; 0) = F,{t) + F^{t), (50) 
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and Paa = Pab = Pba = 0. 

Another interesting quantity that can be calculated, is the long time behavior of this 
functions . The only nontrivial case, is the case (b) (or equivalently (c)). We expect that if 
at t = , we have A and B particles, (one healthy and one diseased individuals), the healthy 
one becomes diseased finally, or we have two B particles finally. In other words, we expect 
that (in case (b)), 

OO 1:2 — 1 

£ PABi^;t^oo\A,B,y;0) ^0, (51) 



and 



X2=y2 xi=yi 

OO X2 — 1 



Y PBBi^;t^^\A,B,y;0)^l. 

X2=y2 xi=yi 

After some calculations , one can show that 



(52) 



OO 2:2-1 



J2 E PAB{^;t\A,B,y;0) = e 



-2t 



2 Y /m(2t) + /o(2t) + 

^y2-yi 



X2=y2 xi=yi 

where In{x) is the n-th order Bessel function of the first kind : 

\n+2k 



{x/2y 



Q k\{n + k)\ 



Now at a; — 00, we have 



In{x) 



therefore, 



00 X2—1 

J2 E PAB{^;t^oo\A,B,y;0) 

X2=y2 xi=yi 



M 

71^' 



(53) 



(54) 



(55) 



(56) 



which goes to zero. M is the number of the In{2tys in the left hand side of (|53|). Using (^81) 
we see that both limits in ( ^Tj) and (|52|) are satisfied. 



4 The class of models 

Now we want to find the class of all two-species integrable models which have the four 
properties introduced in the introduction section. If we look at the preceding sections, 
we notice that all the information about the model are abbreviated in the b matrix (|TB|), 
because this matrix comes from the boundary conditions (^ to (|^), and the latter induce 
our interactions. We also note that the sum of each column of b matrix is one. Now we 
claim that each b matrix which has the following properties: 

1- the non-diagonal elements are one or zero, 

2- the sum of elements in each column is one. 
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represents a model that its interaction(s) can be induced by the master equation (Q) plus 
the boundary condition(s) which can be read from b. 

The reason of the first requirement is that the non-diagonal elements are the sources 
of our reactions, as can be seen from the example solved in section 2.1, and if we want all 
the reactions to occur with equal rate one, the pre-factors of all source terms must be one 
(or zero, if we don't want the corresponding source terms). Note that we should take all 
the rates equal to each other, otherwise for the reactions we are interested in (i.e. those in 
which, particles change their type), the factorization (|n|) will not yield the time independent 
boundary condition(s), which is wrong. 

The reason for the second requirement lies in the conservation of probability. Suppose 
that in the first column of b, for example, the sum of the non-diagonal elements is m and 
the diagonal element is n. So we have m possible interactions each can be a sink for AA. 
Therefore if our configuration is (0AA0), we must have m + 1 sinks, one sink for diffusion 
(0y4A0) (0A0A), and m sinks for reactions: {%AA%) {^a(3^). Now we consider the 
master equation for this process: 

d 



dt 



PAAix^x + l) = Paa{x - l,x + 1) + Paa{x,x) - 2Paa{x,x + 1) 

m' 

= Paa{x -1,x+1) + Y, Papix, X + 1) - (2 - n)PAA{x, x + 1). (57) 



in which we have supposed that the m' elements of the first row of b (besides bn) are 
different from zero and therefore the corresponding probabilities appear in Paa{x,x). Now 
as the number of sinks must bem+1, so2 — n = m + l which yields n + m = 1. Therefore 
the sum of the elements of the first column of b must be one. By the same reasoning, it is 
true for other columns. 

In this way we have 2^^ = 4096 possibilities for matrices b (there are 12 non-diagonal 
elements, each can be one or zero), each plus master equation (|^) can model a reaction- 
diffusion process. But as we have seen in reaction (^, these b's must be consistent with 
the QYBE (|3TD. Therefore the domain of b's is much smaller. So it is sufficient to check 
which of these b's (or more carefully, the S-matrices that are constructed by these b's from 
eq.(p4|)) satisfy (|3T|) . Using a symbolic manipulator (e.g. MAPLE), we found that there are 
42 different b's that satisfy eq. (|3TD . By a closer inspection of these matrices it is observed 
that 14 one of them can be obtained from the others by interchanging A ^ B, so they do 
not represent any new physical interactions. The 42-14=28 b's (interactions) are as follows: 



bi 



/ 1 
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1 














1 





lo 








1 ) 
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It should be mentioned that in above, the reaction processes of each b have been given only 
and the diffusion processes (which exist in all cases) have been suppressed. Also note that 
bi is the pure diffusion process of |]T^, and is nothing but eq. (|1^). 

In all the above cases the probabilities, ^ and A^- are given by (pUD , ( p^Tf ) and (p2D , 
respectively. Obviously, S{z, w) must be calculated from eq. (pif) for each case, and then 
the calculations of section 3 can be repeated for them. 
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